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A  t  .*■  r.r  i  .3^.*  *  r  t’—  :  -rn  i  nj  t  i  on  of  pyrometer  values  in  the  equations 

>  lir':  :  :nvo  ;l  ,,v-\tvrr  is  described  which  uses  a  variation  of  the 
:  i  f  i  ;l  r  r ♦  i  .v-s  c  .ir  vo  -fitting  technique  and  field  data.  Ana- 

I  v  t  |  ■  --  solutions  ure  nof  required,  and  a  digital  computer 

.  ir-  -ai  .  r-;.  *  t  hv  mathematical  solution  to  the  equations  of  the 

"  s  4  ■  I  .  iri,.  examples  in  atmospheric  diffusion  and  a  missile  bel¬ 
li,.,  "  .>*1  Jm/a  how  to  recover  diffusion  parameters  and  aerodynamic 

...  ierfc  ng  field  data  and  the  computer  numerical  solution  to 

t  h**  *  if firr  .j r  i  .  j I  ^ode I  . 
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INTRODUCTION 


The  solution  to  a  mathematical  theory  may  take  a  variety  of  forms  in¬ 
cluding  solutions  which  can  be  described  in  closed  form  or  numerical 
so  I l+ ions  which  requ i re  a  complex  digital  computer  program  if  closed 
form  solutions  cannot  be  derived.  Regardless  of  the  type  of  solution 
obtained,  there  are  usua'ly  a  number  of  unknowns  or  crudely  estimated 
physical  parameter  values.  In  any  case,  some  method  of  finding  the 
set  of  parameter  values  is  needed  such  that  the  model  results  will 
"best”  fit  observation.  One  technique  for  finding  these  ’’best”  values 
is  least-squares  curve  fitting.  Closed  form  solutions  can  be  fitted 
by  one  of  several  methods  to  give  the  best  fit  to  the  data  by  the  least  - 
squares  criterion.  Unlike  the  closed  form  type  of  solution,  complex 
numerical  solutions  are  not  easily  adapted  to  I  east-squares  curve  fit¬ 
ting. 

In  theory  it  should  be  possible  to  invert  the  mcdel  and  use  the  data 
to  evaluate  the  unknown  parameter  values.  One  simple  method  might  be 
to  run  the  computer  program  using  different  input  parameters,  observ¬ 
ing  the  variation  in  the  output  so  as  to  make  a  subjective  estimate 
of  the  best  values.  Repeated  trials  could  narrow  down  the  range  of 
values,  but  the  expenditure  of  individual  effort  and  costly  computer 
time  would  be  undesirable  for  other  than  the  simplest  models. 

This  report  describes  a  technique  which  determines  the  ’’best”  para¬ 
meter  values,  in  a  least-squares  sense,  for  a  computer  simulation  of 
a  complex  physical  process  by  using  the  program,  field  data,  and  a 
variation  of  the  differential  corrections  least-squares  curve-fitting 
method . 


THEORY 

Suppose  that  a  model  of  an  observed  natural  phenomenon  exists  and  that 
to  be  able  to  rely  fully  upon  the  results  that  the  model  can  supply, 
several  parameters  require  numerical  estimates.  One  way  to  determine 
the  best  values  for  these  empirical  constants  is  to  take  data  of  the 
physical  process,  invert  the  model  solution,  and  solve  for  the  unknown 
parameters.  Where  complications  arise  in  the  execution  of  this  method, 
least-squares  curve  fitting  can  be  used  to  recover  the  parameter  values 
which  give  the  ’’best”  estimates  in  a  l east-squares  sense. 

With  data  in  the  form  (xj,  yj),  i  =  1,2,3,  ...N,  one  selects  a  model, 
i.e.,  a  function  F(x;  aj,  02*  •••  a^)  of  x  and  M  undetermined  para¬ 
meters  or,  j  =  1,2,3,  ...M  such  that  F(xj;  3j,  do,  ...  ay|)  will  best 
fit  the  aata  points  yj,  i  =  1,2,3  ...N.  The  ’’best”  set  of  values  for 
the  a.,  j  =  1,2,  ...M,  in  a  I  east-squares  sense  will  be  found  when 


-  me  squares  f  the  differences  in  the  functional  evaluation 
V|  is  j  m ini^ijrr  ,>Ver  all  the  data  points  (Xj,  yj), 

,  .  .  Th.jt  is 


Is  j  'vir'irm~.  "ost  methods  of  least-squares  curve  tiffing  will  speci¬ 
fy  the  f.>r,r  o*  F  whicn  is  ;o  be  used  as  the  model.  The  curve-fitting 
rocnn  i  ■:>  differential  corrections  C  I  — 3l]  requires  no  special  form 
fT  the  "  -be  I  *  unction  F(x.;  a  ^ ,  a0,  ...  a^) .  Restrictions  on  the  pro- 
:ert i es  of  r  are: 


I-  F  must  be  a  s i ng I e-va I ued  function  of  x  and  of  the  M 
parameters  a.. 

2'  r  must  be  differentiable  with  respect  to  the  M  para¬ 
meters  a  . . 

ce cause  differentia  I  corrections  technique  does  not  specify  the  form 
o*  the  ^unction  F  to  be  used,  F  may  assume  the  form  of  a  sequence  cf 
operations  in  a  digital  computer  program.  Most  computer  numerical 
solutions  to  equations  describing  a  physical  process  behave  as  func¬ 
tions  in  that  there  exists  an  independent  variable,  say  x,  such  that 
the  numerical  solution  given  by  the  digital  computer  program  can  be 
regarded  as  a  single-valued  function  of  x,  and  the  dependent  variable(s) 
will  usually  be  d i f ferent i ab I e  in  one  or  more  model  parameters.  It 
should  then  oe  possible  to  use  a  computer  program  as  a  function  F  in 
such  a  way  that  the  methods  of  I east-sauares  curve  fitting  can  be 
applied  tc  F.  This  report  describes  a  method  of  accomplishing  this, 
saving  the  excessive  labor  and  costly  computer  time  involved  in  a 
tr i a l-and-error  approach. 

The  method  of  differential  corrections  and  Its  application  to  the 
I  east-squares  curve  fitting  of  a  computerized  mathematical  model 
soluticr  will  be  discussed  below. 


ver  the  set  of  data  points  (x, ,  y j )  i  =  1,2,3,  ...N,  let  the  differ¬ 
ence  between  the  ^unction  value  F(aj;  aj,  32,  ...  a^)  and  the  data 
pern*  v;  ce  called  the  residual  vj  where 


aM>  " 


(  I  ) 


Mt  cv  least-squares  criterion  will  be  obtained  when  the 

'"in.*-:  '.a  lues  .-j.,  ;  =  1,2,  . .  .M,  give  a  minimum  for 


2 


N 

Z  v 


2 


The  values  for  the  a^s  can  be  found  by  the  solution  of  the  follow! no 
set  of  simultaneous  equations: 


2 

Z  v.  =  0, 
i  =  i 


N 

Z 

i  = 


=  0. 


When  c  is  represented  by  a  digital  computer  program,  the  equations  can¬ 
not  be  solved  directly,  and  an  iterative  scheme  must  be  used. 


Let  Cj  =  a T  +  Aa .  where  the  superscript  o  refers  to  first  approxima¬ 
tions  so  tnat  the  Aajfs  are  the  corrections  needed  to  force  ivj^  to 
a  minimum.  Writing  equation  (I)  with  F(x;  a^  +  Aa j ,  a fi  +  A.:*-,  ... 

+  Aa^)  expanded  in  a  Taylor  series,  we  have 


aF.  aF. 

f-/  O  o  O .  I  I 

V.  =  F(x . ;  a,,  a2,  ...aM)  +  Aa  +  Aa2~ 

3a,  oa, 

[  z 


3F . 

i 


3F . 


+  Aa^ — ~  f  ...  +  A aM- -  -  y.  +  (hiaher  crjer  terms] 

s  r  O  r)  O  I 

daT  raM 

3  M 


Assuming  a  sufficiently  close  estimate  on  the  initial  approximations 
the  higher  order  terms  will  be  neglected.  Defining  R.  as  the  residua! 
Rj  =  Yj  -  F?  where  F<?  =  F(x.;  a®,  ...  aj^)  equation  { :>  c  sr>  bo 

written 


v . 

i 


c*F 


3F. 


3F. 


Aa  . -  +  Aa-v -  +  Aa  , -  + 

L  o  2,o 

'?a| 


I' 


The  system  of  equations  (2)  which  are  to  bo  solved  r/wr  * 
equat i on 


3F. 

9F. 

9F. 

9F. 

I(Aa;-C 

t  Aa  , — -  +  ... 

+  Aa.; - 7  - 

M 

R.)  — -  -* 

i  '3a, 

o 

1  9a° 

J 

Letting 

9F.  9F. 

V  ■  f 

i  i 

.  0  7o 

9a  .  oa 

J  * 

and 

9F. 

B  .  =  I 

R.  — -  , 

i 

1  9a? 

J 

we  have 

the  algebraic  system  of  s 

i  mu  1 taneous 

equations 

matr i x 

form 

/A 1  1  A 1 2 

A  1 3  •••  AIm\ 

f\  ■ 

M 

A2 1  A22 

A  A 

23  2M 

Aa2 1  ' 

(5) 


Ml 


M2 


M3 


MM, 


fW 


\ 


\V 


Since  Aj#k  =  .  the  matrix  system  above  is  symmetric  allowing  only 

"■fie  evaluation  6^  the  upper  triangle  and  making  use  of  the  square-root 
^ethcd  for  solving  symmetric  systems  of  simultaneous  equations  El, 4,5]. 


Having  obtained  the  corrections  Aa:,  j  =  1,2, 3,  .  .  .M,  they  either  satis¬ 
fy  the  least-squares  requirement  when  added  to  a?,  giving  the  a 


A a  : ,  or  the  newly  formed  a 


^  u,,  1  '  '-J 

■ 1  s  give  a  new  set  of  first  approximat 


=  a  3 


-  -  3  '  '  ”  -  - - -  ~  ’  1  ’  1  -  '  rr  - -  ’ 

say  Jnd  a  new  set  of  corrections  Ab j  is  evaluated.  Convergence 
f or  +he  method  may  be  said  to  have  been  reached  when  no  significant 
chance  in  the  set  of  model  parameters  a j ,  j  =  1,2,  ...M  is  exhibited 


ions 


tn at  is,  when 


*  results  1 rom  cycle  to  cycle. 
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To  apply  the  method  of  differential  corrections  curve  fitting  to  a 
computerized  numerical  solution,  one  must  use  a  numerical  approximati 
for  each  of  the  required  partial  derivatives 


3F. 

l 


as  defined  above.  Use  is  made  of  a  seven-point  numerical  formula 
described  below 

Given  discrete  points  (aj,  F. )  of  a  tabulated  function  Fj  =  F ( a , ) 
for  i  =  -J,  -(J-l),  ...  -I, 0,1, 2,  ...  J  for  equispaced  points  a.  with 
a^+j  -  a |  =  h,  the  first  derivative 


evaluated  at  aQ  is  given  as  follows: 

*  aO  -  f-i  -  Xfr  V,*  *  w  <0  -  ‘V’j 

with  6 F |  =  Fj  +  |^2  "  ^|-|  ?•  This  gives  the  formula 

F1  =  [F.  -  F  .  -  l/6(F0  -  2F.  +  2F  ,  -  T 

o  I  - 1  /  I  - 1  -z 

*  IiJ,F3  -  F2  -3(F2  -  Fl>  *  3<FI  -  Fo>  -  ro  *  F-l> 

-  3B<FI  -  Fo  -3<Fo  "  F-l>  *  3<F-I  -  F-2'  -  F-2  * 

Use  of  this  algorithm  gives  the  required  partial  derivatives  for  the 
differential  corrections  method  described  above.  However,  one  may 
improve  on  the  evaluation  of  partial  derivatives  in  the  event  that 
one  or  more  of  the  model  parameters  a-  are  used  in  a  sub  function  f  of 
the  numerical  routine.  When  these  individual  functions  f  (i.e.,  one 
of  the  governing  equations  in  the  me»themat  i  ca  I  formulation  of  the 
o^oblem)  can  be  defined,  it  is  recommended  that  The  chain  rjlo  for  ca 
culating  derivatives  be  used.  Letting  F  be  the  main  program  solute  o 
and  f  an  intermediate  solution,  a  subprogram,  or  the  like,  thor, 


a  a  .  9  f  9a, 

J  J 

This  rule  can  provide  better  convergence  ir  the  r « >u 4  i  r,.  • 


EXAMPLES 


Atmospheric  Diffusion 

An  atmospheric  diffusion  model  for  the  prediction  of  contaminant  con¬ 
centrations  in  the  atmosphere  will  be  treated  with  the  curve-f i tt i ng 
technique  previously  described.  The  concentration  C  is  a  function  of 
downwind  travel  distance  x  and  several  atmospheric  parameters,  together 
with  empirically  determined  cons+ants.  The  mathematical  theory  used 
is  based  upon  the  application  of  the  similarity  theory  [7,83  of  the 
surface  boundary  layer  as  related  to  diffusion  in  the  atmosphere  [9]. 

The  concent  rat i on  function  C(x;  a^,  a2>  a^)  is  given  by 


2  2 

L  u*  ;  <f(r,>  -  f(c  )) 
*  o 

where 


dx  u* 

(fU)  -  f(co))' 
rr  =  bu„  0<  O, 


0(5)  -  (.  -  f 


-I 


1/4 


f(p)  =  ln(cc/pQ)  +  6(C-C  > 
where  the  sunbols  are  defined  as  follows: 

5  =  source  strength 
L  =  atmospheric  parameter 
u*~  atmospheric  parameter 
k  =  K-1rm3n?s  constant 
i  r  i  ca I  constant 
:al  constant 

vj  travel  distance  of  substance 


z  =  vertical  height  coordinate 
t  =  time. 

These  representative  equations  for  modeling  atmospheric  diffusion  urine 
similarity  theory  have  been  solved  numerically  using  a  digital  compu¬ 
ter.  The  curve-f i tt ing  technique  previously  described  will  be  tested 
by  evaluating  the  empirical  constants  c  and  8,  and  also  by  inserting 
into  the  0  function  two  parameters  to  be  fitted,  a j  and  given  by 

0U>  -  «.,  -  jifi  ')  “2. 

Original  speci f i cat  ions  for  the  a’s  are  =  I  and  =  1/4. 

For  this  example,  a  test  case  was  developed  by  specifying  c,  3,  otj, 
and  c*2  to  generate  a  sample  data  deck.  Then,  approximations  to  c, 

8,  oi[,  and  0L2  are  given,  and  by  using  the  data  the  original  parameters 
will  be  recovered . 

First  only  c  and  8  are  evaluated  as  shown  in  Table  l-A.  Then  c,  8, 
a | ,  and  a 2  are  evaluated  together  as  given  in  Table  l-B,  The  results 
of  the  test  case  show  how  model  parameters  for  atmospheric  diffusion 
theories  can  be  evaluated  using  field  data  and  the  curve-f i tt i ng  tech¬ 
nique  described  herein. 


A  Ballistic  Missile 


The  curve-fitting  technique  was  applied  to  the  problem  of  determining 
the  drag  coefficient  of  a  ballistic  missile  using  radar  position  data 
of  the  missile  during  flight  [l].  To  use  the  data  in  a  curve  fit  of 
the  type  described,  an  accurate  realistic  model  of  the  trajectory 
during  all  phases  of  the  missile’s  performance  is  required. 

The  ballistic  model  selected  is  a  f i ve-degree-of-f reedom  ballistic 
missile  trajectory  and  impact  prediction  digital  computer  program, 

TRAJ  [10],  used  for  missile  flight  safety  at  White  Sands  Missile  Range, 
and  also  utilized  in  studies  of  meteoro log i ca I  effects  on  a  ballistic 
missile  by  the  Atmospheric  Sciences  Laboratory.  Only  the  three  posi¬ 
tion  coordinates  relative  to  the  launcher  will  be  used  in  the  test 
case.  Hence,  X,  Y,  Z  are  the  dependent  variables*,  and  Mach  number,  m, 

? - 

Differential  corrections  method  applieu  to  a  vecf or-va I ued  function 
is  described  in  Appendix  A. 


Final  residual  R  =  .0973339 


C-97??^ 


I  *  ;  .  ^c-c  to  be  the  only  independent  variable  of  the  model, 

f  vp.  jru  m.jny  input  parameters  to  TRAJ,  only  the  aerodynamic 

’■**!  !■•■  r:  -f  the  drag  curve  will  be  considered  as  influences  on  the 

,  -  i  It-'  .  :r  j  ■o^t'.-ry. 

■*  <v . i  i  g^r’i^ratec  by  using  a  known  d^ao  Curve  is  input  and  put- 

■,  />  ..^tput  on  data  cards.  The  input  drag  curve  was  con- 
r-  ’--J  H:  :ivt‘  the  character  i  st  i  c  shape  and  values  of  an  actual  drag 
'  .•  ,  :l:r-  j.:h  most  drag  curves  are  provided  in  the  form  of  a  table. 

.* ’  he ! no  .ould  "piecewise”  fit  for  a  drag  function  if  no  single 
sroS:'.  i .  ’  can  :e  used.  Hence  the  drag  curve  used  is 

A 

t  x  j  i  5 .  ,  ( m+  |  )  ** 

o  --  •  i  +  a.,  x  +  x ( m  +  .8)  )  Lxp(-  — — - ) 

Cl  do  10 

i*  it i ally  values  for  ct|  ,  X2>  and  a-^  are  .25,  .5,  and  .75,  respec- 
/ely.  Results  of  the  iteration  process  are  given  in  Table  II.  The 
rsr  arprox imat ic ns  were  specified  and  the  original  values  of  the 
rjmeters  were  obtained,  to  a  certain  degree  of  accuracy. 

SUMMARY 

is  report  describes  a  method  of  curve  fitting  for  a  broad  class  of 
notions  not  representab I e  by  simple  terms  in  a  closed-form  expression, 
est  ^unctions  can  assume  the  identity  of  complex  sequences  of  opera- 
•ns  in  a  digital  computer  program  and  yet  have  suitable  properties 
r  o^rve  fitting.  This  technique  will  afford  considerably  more  f I  ex¬ 
ility  to  the  engineer  or  scientist  who  wishes  to  curve  fit  a  function 
*ield  data  without  overly  simplifying  the  function  to  be  used.  The 
vest i gator  can  now  build  a  complex,  sophisticated  mathematical  model 
u  physic.: I  process,  solve  the  equations  numerically  on  a  digital 
mooter  and  use  field  data  and  the  method  described  here  to  find  un- 
own  par reters  cf  the  physical  system. 
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APPENDIX  A 


Suppose  data  have  been  observed  in  which  there  ex  i r  •  s  rv-re  than  r.t  -  ;«- 
pendent  variable  for  a  common  independent  variable.  i  hoc  the  r  j . :  i . :  i  |  : 
be  of  the  form  (x*,  Y.)  where  '2 •  is  a  vector  of  L  components.  Lei 
FCxj;  aj,  62 ,  ...  a^)  be  a  vector-val  ued  function  of  L  .  ■•re  r-er:  t  a  whir/ 
is  to  approximate  the  Y;  in  a  least-squares  sense.  Iso  pr.blor  is  h 
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the  solution  of  which  is  identical  to  the  solution  given  in  the  text, 
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